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1.   INTRODUCTION 

It  has  long  been  appreciated  that  the  simple  birth-and-death  Markov 
process  often  provides  an  adequate  initial  model  for  the  behavior  of 
service  systems,  populations,  epidemics,  and  many  other  stochastic 
systems;  see  Feller  [3]  for  an  early  classic  account.  Refinements, 
particularly  in  the  modeling  of  service  systems,  have  typically  involved 
the  replacement  of  exponential  service  times,  also  assumed  independent,  by 
independent  random  variables  of  "general"  or  arbitrary  distribution, 
replacement  of  independent  exponential  inter-arrival  times  by  independent 
"general"  random  variables,  or  both.  Isolated  examples  in  which  arrival 
process  parameters  are  allowed  to  change  deterministically  in  time  have 
also  been  studied. 

Randomly  appearing  fluctuations  in  system  environment,  associated 
with  weather  or  other  change  in  physical  surroundings,  personnel  changes, 
alteration  of  system  usage  intensity,  etc.,  are  likely  to  be  reflected 
in  changes  in  observed  failure  and  repair  rates  in  a  system  reliability 
context. 

1.1  Examples  :  In  order  to  illustrate  the  ideas  summarized  above,  we 
introduce  two  specific  birth  and  death  models  in  random  environments. 
Numerical  results  and  discussion  will  be  found  in  the  last  section. 

Repairman  Model .  Suppose  any  of  m  machines  that  are  in  use  fail  independent- 
ly in  Markovian  fashion  at  rate  {A.,  t  >  0},  and,  if  failed,  experience 
repair  at  rate  {p.,  t  >  0}.  In  turn,  the  rates  {*.}  and  {p.}  are  then- 
selves  finite-state  Markov  processes  independent  of  the  state  of  the 
machines,  so  that  if  n  is  the  number  of  machines  on  repair,  and  j  identifies 
the  environmental  state,  then  (n,j)  is  the  overall  state  variable  of  the 


system;  conditional  upon  j,  n  changes  by  one  unit  at  a  time  in  typical 
birth-death  fashion  with  parameters  depending  on  j  :  the  probability  that 
an  operating  machine  goes  down  in  (t,t+dt)  is  x,(m-n)dt  +  o(dt),  while  tne 
probability  that  a  machine  on  repair  becomes  available  is  y  .  min(R,n)dt  +  o(dt)  , 

being  the  number  of  repairmen.  Such  a  setup  describes  groups  of 
redundant  equipments  that  all  experience  common  environmental  intensities 
simultaneously;  the  environmental  changes  are  reflected  in  the  numerical 
values  of  the  failure  and  repair  rates  that  prevail  at  any  time  point. 

The  model  described  also  may  be  used  to  represent  the  behavior  of  a  net- 
work of  m  timesharing  computer  terminals  that  independently  send  messages  or 
programs  to  a  common  central  computing  facility.  Suppose  the  rate  of  message 
transmission  changes  with  external  activity,  i.e.  responds  to  an  occasional 
period  of  unusual  activity,  perhaps  a  crisis  situation.  In  this  case  the 
"constant"  demand  rate  A  (assumed  equal  for  all  terminals)  switches  almost 
instantaneously  to  a  higher  value,  switching  back  to  normal  after  the  crisis 
elapses,  but  doing  so  repeatedly.  Left  alone,  the  central  processor  may  well 
continue  processing  at  the  original  rate,  allowing  congestion  to  simply  build 
up,  and  eventually  drop  again  when  the  demand  lapses.  On  the  other  hand, 
remedial  action  may  be  taken.  These  phenomena  suggest  interesting  and  realis- 
tic questions  concerning  stochastic  control,  but  these  will  not  be 
considered  in  this  paper. 

Mass  Search  Model .  Let  a  group  of  m  predators  attempt  to  round  up  and  capture 
a  finite  group  of  p  prey.  Predators  search  independently  and  with  equal  inten- 
sity and  effectiveness,  but  the  detection  rate  is  allowed  to  depend  upon  exter- 
nal environmental  conditions  that  cause  relatively  long-term  changes  in,  say, 
visibility.  If  a  prey  is  detected,  it  is  followed  until  lost  by  the  predator, 
after  which  moment  it  is  susceptible  again  to  search,  detection,  and  active 
surveillance.  While  predators  and  prey  move  independently,  as  soon  as  a  pre- 
dator begins  following  a  particular  prey  the  latter  is  removed  from  circula- 
tion and  the  remaining  unattached  predators  continue  search  for  the  free  prey. 


Of  interest  is  the  long-run  or  stationary  distribution  of  the  number  of 
prey  under  simultaneous  surveillance,  and  also  the  first-passage  time  until 
a  large  fraction-  perhaps  all-  of  the  prey  are  simultaneously  under  the  eye 
of  the  predators.  Here  is  a  plausible  model  :  the  state  of  the  system  is 
(n,j),  so  that  if  the  number  of  prey  under  surveillance,  is  n,  and  the 
state  of  the  environment  is  j,  then  the  probability  that  a  free  prey  is 
detected  is  A  .(p-n)(m-n)dt  +  o(dt),  while  the  probability  that  a  prey  under 
surveillance  is  lost  is  y.ndt  +  o(dt),  where  n  <  min(p,m)  and  \.   and  u 

J  J  J 

change  in  accordance  with  independent  finite-state  Markov  processes,  as 
before. 

1.2  The  analytic  structure  :  We  shall  consider  Markov  processes 
{(Xt,Yt),  t  >  0},  on  the  state  space  {(n,j),  0  <  n  <  N,  1  <  j  <  Kn),  with  a 
block-tridiagonal  infinitesimal  generator  Q  : 


A<°>   A(°>  0     0 

nW   a*1)  a*1)  0 

0     M<2)  A<2>  A^ 


,(N) 


M(N-1)  A(N-1)  A(N-1) 


,(N) 


(1.1) 


Here  A^  ;,  A^   ,  ...,  A1    are  square  matrices,  respectively  of  order 

Their  diagonal  elements  are  strictly  negative,  all  other 
elements  are  non-negative.  The  matrices  a'  ,  0  <n<N-l,  and  M^, 

=  N,  are  rectangular,  with  appropriate  dimensions  ;  their  entries 
are  non-negative.  The  rowsums  of  Q  are  equal  to  zero,  therefore  we 
have  that 

A(°)e+A(°)e  =  0  . 

M(l)e  +  A(l)e  +  A(i)e  =  0,     1  <  i  <  N-l, 

M^e  +  A<N>e  -  0, 

where  e  denotes  a  column  vector  with  unit  elements.  The  variable  Y 
is  to  be  interpreted  as  the  state  of  the  environment,  and  X   as  the 
state  of  the  birth  and  death  process,  at  time  t  . 

We  assume,  furthermore,  that  the  Markov  process  Q  is  irreducible,  and  we 
denote  by  level  n  the  set  {(n,j),  1  <  j  <  K  },  of  states  corresponding  to 
the  common  value  n  for  the  first  index.  The  structure  (1.1)  of  Q  permits 
the  Markov  process  to  move  up  or  down  by  only  one  level  at  a  time.  It  is 
in  this  respect  analogous  to  the  classical  birth-and-death  process,  see 
Feller  [3]  or  Karlin  and  Taylor  [5]. 

In  the  examples  cited  above,  either  X  or  Y  changes  each  time  the 
Markov  process  undergoes  a  change  of  state.  In  fact,  this  restriction  is 
not  part  of  the  model,  and  we  allow  X  and  Y  to  change  simultaneously. 

1.3  Existing  litirature:  Infinite  birth-and-death  models  in  random 
environment  have  been  studied  for  some  time  already;  early  results  are 
found  in  Eisen  and  Tainiter  [2],  Purdue  [11]  and  Yechiali  [14].  More 
recently,  Neuts  [9],  Chapter  6,  has  systematically  examined  a  class  of 
problems  in  which  N  =  »,  and  Q  has  a  special  repetitive  structure  : 

A<»  .   A<2>=  .... 

a'1).  a<2>=..., 

M<2)  »  M<3)=  .... 


This  structure  leads  to  matrix-geometric  stationary  probability  vectors  and 
efficient  algorithmic  procedures. 

Finite  models  have  been  examined  by  Torrez  [12],  who  suggests  using 
numerical  procedures  designed  to  solve  eigenvectors  for  band-matrices. 
Hajek  has  considered  in  [4],  Section  5,  a  finite  model  with  repetitive 
structure  : 


A<D 

=A<2>=. 

,.=A<N-2>. 

A*1) 

=  A<2>  =  .. 

..■a'""1). 

Ml2> 

=M<3>=. 

..  .  mC-1), 

Kl 

-h       •- 

..    =  KN_j   =   K. 

Hajek  determines  the  stationary  probability  vector  in  terms  of  two  matrices 
R  and  R,  of  order  K,  which  have  to  be  iteratively  computed.  Finally,  Keilson 
et  al .  [6]  have  considered  finite  models  with  the  structure  (1.1)  and  equal 
K  's.  They  analyse  the  (Laplace  transform  of)  first  passage  times  distributions, 
from  which  they  obtain  equations  for  moments,  and  for  the  stationary  distribu- 
tion. We  cannot  in  this  short  space  describe  in  detail  the  differences  between 
our  results  and  those  in  [4,  6,  12];  but  we  shall  give  some  discussion  at 
the  end  of  Section  5. 

This  paper  presents  an  efficient  computational  approach  to  the  analysis 
of  birth-and-death  models  in  a  Markovian  environment.  The  emphasis  is  upon 
obtaining  numerical  properties  of  both  stationary  distributions  (in  the  next 
section),  and  first-passage  time  (in  Sections  3  and  4).  The  computational 
algorithms  are  discussed,  and  numerical  examples  are  given,  in  the  last  two 
sections. 

2.  THE  STATIONARY  DISTRIBUTION 

In  order  to  determine  the  stationary  probability  distribution,  and  moments  of 
first  passage  times,  it  is  useful  to  think  of  the  Markov  process  as  evolving  in 
a  certain  manner. 


For  0  <  n  <  N-l,  define  S   to  be  the  restriction  of  the  original  process  Q, 
observed  during  those  intervals  of  time  spent  at  level  n,  before  the  original 
process  enters  level  n+1  for  the  first  time,  the  state  space  of  Sn  is 
{(n,j),  1  <  j  <  K  }.  Clearly  all  S  ,  0  <  n  <  N-l  are  transient  Markov  proces- 

The  process  SN  is  the  restriction  of  the  process  Q  to  the  states 
{(N,j) ,  1  <  j  <  K.,};  it  is  an  ergodic  Markov  process.  We  denote  by  C  the 
infinitesimal  generator  of  the  process  S  ,  0  <  n  <  N. 
In  order  to  determine  the  matrices  C  ,  we  need  the  following  result. 

Lemma  1 

Consider  a  Markov  process  on  the  state  space  {1,2, . . . ,r,r+l,r+2, . . . ,r+s}, 

with  infinitesimal  generator  q  : 


% 


0 


(2.1) 


where  A  is  a  square  matrix  of  order  r,  A  is  a  rectangular  r  by  s 
matrix.  The  states  r+1  to  r+s  are  all  absorbing. 

are  all  transient  if  and  only  if  the  matrix  A  is  non- 
singular. 

b.  The  (i,j)th  entry  of  (-A-1),  for  1  <  i ,  j  <  r;  is  the  expected  amount  of 

in  the  transient  state  j,  starting  from  the  transient  state  i, 
before  absorption  in  any  of  the  absorbing  states. 

c.  The  (i,k)th  entry  of  (-A"  a),  for  1  <  1  <  r,  r+1  <  k<r+s,  is  the 
probability  that,  starting  from  the  transient  state  i,  absorption  occurs 
in  the  state  k. 


Proof.  The  first  assertion  is  proved  in  .Neats  [9],  Lemma  2.2.1. 
The  matrix  A"  has  all  nonpositive  entries.  The  second  assertion  is 
a  consequence  of  Neuts  and  Meier  [10],  Corollary  2.  To  determine  the 
probability  PH  ,,  of  being  eventually  absorbed  in  state  k,  r+1  <  k  <  r+s 
starting  from  state  i,  1  <  i  <  r,  we  study  the  Markov  chain  embedded  at 
instants  immediately  following  a  transition  in  the  Markov  process. 
The  corresponding  transition  probability  matrix  P  is  given  by 


P  = 


I  -A_1A 


0 


-A^A 


I 


where  the  matrix  A  is  diagonal,  with  diagonal  entries  equal  to  those  of  A. 
It  results  from  Kemeny  and  Snell  [7],  Theorem  3.3.7,  p.  52,  that 

P  =  [I  -  (I  -A_1A)]_1  (-A_1A)  =  (-A_1)A; 
which  completes  the  proof  of  Lemma  1.  i 


We  may  now  determine  the  matrices  Cn. 
Lemma  2 


C  -  A<°> 


C  =AW  +M(n)(-C-\)A(n_1).  Kn< 
n  n_  x 


(2.2) 
(2.3) 


Proof.  The  equation  (2.2)  is  obvious  starting  from  level  0,  the  process  SQ 
terminates  as  soon  as  the  process  Q  enters  level  1.  Meanwhile,  the  transi- 
tions are  governed  by  A(0).  Since  the  process  Q  is  irreducible,  the 
process  SQ  contains  only  transient  states,  and  the  matrix  CQ  is  non- 
singular,  by  Lemma  1.  Also,  (-CQ)  >0. 


We  may  now  determine  the  generator  C1  of  the  process  S^.  Let  Z^x), 

t  >  0,  be  defined  as  follows.  Z^x)  =  i  if  the  process  S1  is  in  state 

(l,i)  at  time  x.  Assume  that  Z^x)  =  i,  and  Z^x+  dx)  =  j  *  i. 

Since  j  *  i ,  a  transition  must  have  occurred  in  the  process  Q.  Z-^x  +  dx)  =  j 

if  and  only  if  one  of  the  following  events  have  occured. 

a.  The  transition  is  from  (1,1)  to  (l,j),  this  happens  with  probability 
aJ1!  dx  +  o(dx). 

b.  The  transition  is  from  (l,i)  to  (0,k),  for  some  k,  a_nd  the  process  Q 
returns  to  (l,j),  after  spending  an  unspecified  amount  of  time  at  level  0, 
before  visiting  any  other  state  at  level  1.  This  happens  with  probability 
Mi  /  dx  .  y.    •  +  o  (dx),  where  y.     ■     is  the  probability  of  moving  from 

1  ,  K         K ,  J  K  ,J 

(0,k)  to  (l,j)  before  visiting  any  other  state  at  level  1. 

If  we  set  in  (2.1)  A  =  A^  '  and  A  =  A^  ',  it  results  from  Lemma  1  that 

It  is  now  clear  that 

=  JlZ^x)  =  i]  =  a{|]  dx  +  1°  m||)  (-A(0)"1A(0))kj  dx 

+  o  (dx),  1  <  1*  j  <  Kp 

and  similarly,  that 

PEZjd+dx)   =   ilZ^x)   =   i]   =   (1+AJJj   dx)   +     1°  m(}>    ("A^'V0*)^.   dx 

+  o   (dx),  1  <  i  <  K1  . 


Thus.     C1=a(1)+m(1)(-,(0))-VO), 

.  ^h  m(1»(-c3V0), 

which  proves  (2.3)  for  n=l. 

Assume  that  (2.3)  holds  for  n,  we  prove  now  that  it  holds  for  n+1.  Let 
Z  ,  (t),  t  >  0,  be  equal  to  i  if  the  process  S  .  is  in  state  (n-t-l.i) 
at  time  t  . 

If  Z  .  (t)  =  i,  then  Zn+1  (x+dx)  =  j*i  if  and  only  if  one  of  the  fol- 
lowing events  occurs. 

a.  There  is  a  transition  from  (n+l,i)  to  (n+1, j ) ,  with  probability 

A^t^dT  +  o(dx). 
■  »J 

b.  There  is  a  transition  from  (n+l,i)  to  (n,k)  for  some  k,  with  probability 
Mi  .  'dx  +  o(dx),  and  the  process  Q  returns  to  (n+1, j ) ,  after  spending  an 

1  ,  K 

unspecified  amount  of  time  at  levels  0,1 n,  before  visiting  any  other 

states  at  level  n+1,  with  probability  y>    \. 


Since  S   records  the  visits  made  by  Q  at  level  n  before  reaching  level 
n+1,  we  have  by  Lemma  1  that 


It  is  now  easy  to  prove  that  (2.3)  holds  for  n+1,  which  completes  the  proof 
of  the  lemma . 


u. 


We  shall  now  determine  the  stationary  probability  vector  P,  i.e.  the 
solution  of  the  system  P  ()  =  0,  P  e  =  1.  We  partition  that  vector  as 

L  =   (fn' — 1'  "'•'  -^P '  wnere  tne  subvectors  P   have  K   elements  and 
correspond  to  the  states  of  level  n,  0  <  n  <  N. 

Theorem  1. 

The  vectors  P  ,  0  <  n  <  N,  are  determined  by  the  equations 

PNCN=0>  (2.5) 

L)(-C^),  for     0<n  <N-1,  (2.6) 


N 

I    P    e    =  1.  (2.7) 

n=0         ~ 


if.  From  the  structure  (2.1)  of  Q,  it  results  that  the  system 
j^  Q  =  P_  may  be  decomposed  into 

Pq  A(0)+  Pj  M(1)  =  0,  (2.8) 

•)+  P  A(n)+  p    M(n+1)  =  0,  Kn<N-l.    (2.9) 
— n+1        —  ' 


•^A^U.  (2.10) 


The  equations  (2.8)  to  (2.10)  are  matrix  equivalents  of  the  familiar 
stationary  equations  for  birth-and-death  processes.  They  may  also  be 
viewed  as  equations  of  probability  balance. 


11. 


From  (2.9),  Pq  =  -Pj  M(1) ( A(0) )_1=  P^  M(1)(-C~V  Then,  by  recurrence, 
using  Lemma  2.  the  vectors  P^  ,  1  <  n  <  N,  satisfy  the  equations  (2.5)  and 
(2.6).  Since  the  matrix  C^  is  the  infinitesimal  generator  of  a  finite 
irreducible  Markov  process,  Equation  (2.5)  has  a  unique  solution,  up  to  a 
multiplicative  constant,  and  that  constant  is  determined  by  (2.7). 

a 

This  result  suggests  an  algorithm  to  compute  the  stationary  probability 

vector. 

Algorithm  A. 

Al.  Determine  recursively  the  matrices  C  ,  0  <  n  <  N. 

A2.  Solve  the  system  tt^  C^  =  0,  n^  e  =  1. 

A3.  Compute  recursively  the  vectors  P^  ,  n  =  N-1,  ....  0,  using  j^ 

instead  of  P^  . 
A4.  Re-normalize  the  vector  P^  so  obtained. 

A  complete  analysis  of  this  algorithm  is  deferred.  We  shall  make  some 
comments  on  it  at  the  end  of  this  section.  Examples  of  numerical  appli- 
cations to  certain  specific  models  are  presented  at  the  end  of  the  paper. 

3.   FIRST  PASSAGE  TIME  TO  HIGHER  LEVELS 

We  denote  by  T   the  first  passage  time  from  level  n-1  to  level  n, 

and  by  x    the  first  passage  time  from  level  n-1  to  level  m  >  n  : 
J       n,m  r    a 


Tn  m  =  inf  {t  >  °  :  Xt  =  m'  X0  =  n_1}  »  l   <  n  <  m  <  N' 


n  "  Tn,n  • 


We  define,   for     x>0,    l<n<N,    l<i<K     , ,   1  <  j  <  K   . 

n-1  J         n 


gJ;J(x)   -  P   [Tn<x,   YT+0  =  j    |    X(J-n-l,   YQ  =   1] 

n 


12. 


the  state  of  the  system  at  time  h,  a  simple  probabilistic 
argument  shows  that 

g(")(x+h)  .  (HA(n:l>h)g(n)W  ♦  I  A^'h  g<°j(x)  ♦  ^h 

•  »J      k*  i 

K         K  (3  1) 

+ Sfa1)h  £i9&1)m  9^(x) + o(h)- 


denotes  the  Stieltjes  convolution.  Substracting 

es  of  equation  (3.1),  dividing  by  h  and  letting 

:irst  order  differential  equation  for  g  n  (x). 

•  »J 

'.(C)  the  Laplace-Stieltijes  transform  of  g^Ux), 
i  ,j  i  >j 

■  gjnl(x)]   dx, and  by  G(n^(C)   the  matrix  with  entries 

sg<o. 

We  conclude  from   (3.1)     that 

C  G(n)(c)=A(n-1)G(n)(C)tA(n-1)   +  M^G^U^^c), 

for  2  <  n  <  N, 
and  similarly, 

tS^'tO-AWs'^lO+iW.  (3.3) 


(3.2) 


These  equations  may  be  written  as 

G(1)(C)   =  DQ(0   A(0),  (3.4) 

Dn-1(0   A(n_1),               for  2<n<N,                                 (3.5) 

where           DQ(C)   ==   (d   -  A(0))_1,  (3.6) 

-  A(n)-  M(n)6(n)(0)_1.  1  <n  <N"1"                       (3-7) 


13. 

Let  gft^x)  =P[Vm<x,  Y^  =  j  |  XQ  =  n-1,  YQ  =  i]  , 

and  let  G(n!m)(£)  denote  the  Laplace-Stieltjes  transform  of  q(n!m^(x). 
We  readily  obtain  that 

G(n.m)(5)  .  s(n.m-l)(0  fi(m)(0i  (3>8) 

m  (k) 
=  M  G^K)(0,  (3.9) 

k=n 

for  1  <  n  <  m  <  N,  where  we  define  G^n'n"^(c)  =  I  for  all  n. 

We  easily  prove  from  Lemma  2  and  Equations  (3.4)  to  (3.7)  that 

Dn  (°)  =  "  Cn1  '  for  0  <  n  <  N-1,  (3.10) 


hence 


:(n) 


(0)  =  -  C"^  A(n_1),     forl<n<N,  (3.11) 


:(n) 


which  is  merely  another  representation  of  Equation  (2.4)  ,  since  (G^  (0)).  . 
is  the  probability  that,  starting  from  (n-l,i),  the  process  Q  visits 
(n,j)  before  visiting  any  other  state  at  level  n. 

Let   u(n'm>  =  -  -  G(n«m)(Olf  n  • 
We  have  that 


u(n,m) 


T"E[Tn.-'YT   +0=J  I  X0=n-1,  YQ=i]  , 

n,m 

and    ^-ElTj^n-l,  Y0-1]  , 

where  u(n'm>  =  U(n'm)e. 

We  define  U(n>  =  U(n'n>  =  -  -   G<n>(Ol,  n  ,  and  u(n>  =  U(n>< 

H        5  u 


14. 


Theorem  2 


xpected  values  of  the  first  passage  time  to  higher  levels  satisfy 
following  recurrence  relations. 

u^  -u*1'1)  —  C^e,  (3.12) 

»  _  u(n,n)  s  _  C-l  (e+M^"1^^'15),   for  2  <  n  <  N,  (3.13) 
y.    _        n-lv—     — 

1)  +  G(n'm-1)(0)  u(m),     forl<n<m<N.  (3.14) 


n  may  be  proved  by  differentiating  equations  (3.6)  and 

and  evaluating  the  derivatives  -r-  D  (c)  at  E,  -   0 
;  from  first  probabilistic  principles.  First,  equation 
lows  from  Lemma  1  (b)  and  (2.2).  Then  let  V  be  the  elapsed 
time  until  the  process  enters  a  different  level;  that  is, 

X.  #  Xq}  .  By  the  strong  Markov  property,  Lemma  1 ,  and 
Lemma  2 

E^Tn+llX0=n  •  Y0=  i] 

:Vl  ^  Xv  =  n  +  1|XQ  =  n  ,  YQ  =  i] 


Xy  =  n  -  1|XQ  =  n  ,  YQ  =  i] 


Kn-1 
e(i)  +  I     [(-A(n))_1  M(n)]ik  E[Tn|XQ  =  n  -  1  ,  YQ  =  k] 

K—  I 


)DE(-C-11)  A^'^^jECT^IXq  =  n  ,  YQ=  j]       (3.15) 


Hence, 


[I  -  (-A^r1  M(n)(-C-]1)A(n-1)]u^n+1) 


15. 


=  (-A(n))-]  e  +  (-A^)"1  M(n)  u(n) 


(3.16) 


Multiplying  both  sides  of  (3.16)  by  A^  results  in  the  equation 


[A(n)  +M{n)(-C-\)An-h   u(n+1) 
n- 1      — 


(3.17) 


+  M(n)  u(n)j 


Equation   (3.13)   now  follows  from  (2.3)   and  (3.17). 
Since 


EFVm-l   *  TJXj,  -  n  -  1 .  YQ  -  1 D 


=  U(".m-D(i)   +G(n'm-1)(0)   u(m)    (i) 


equation  (3.14)  also  follows. 


Any  number  of  moments  of  the  first  passage  times  may  be  similarly  obtained, 
We  merely  state  the  following  result  for  second  moments,  without  proof. 


10 


Le 


t     vin'm)  =  E^n,m!  X0  =  n-1,  Y0  =  i], 


v(n)  =  v(n>n)5 
for  1  <  n  <  m  <  N.  We  have 

v(i)  -  2  r-  c"1^1) 

±        ~  d   (  L0  >±       '  (3.13) 

„(n)   9  /  r"l  \  m   M(n-l)n(n-l)\  (n) 

v     2  ("  cn-r  (T  +  M    u    )±  (3.19) 

+  ("  Cn".\)  H1""1^"1),   for  2<n<N, 
,(n,m)  =  ^(n,m-l)  +  2  ^n.m-iym)  +  s(n,m-l)(0)  ^(m),         (3>2()) 

for  1  <  n  <m  <  N. 


s.  Theorems  1  and  2  show  how  the  matrices  C  ,  0  <  n  <  N,  determined 

2,  play  a  central  role  in  the  determination  of  both  the  stationary 

ility  distribution,  and  the  moments  of  first  passage  times  to  higher 
levels . 


16. a 


Example:  Consider  the  repairman  model  of  Section  (1.1).  Assume  there 
are  just  two  environment  states,  denote  by  j  =  1,  2  .  Let  the  transition 
rate  from  environment  state  2  to  1  be  a  ,  and  from  environment  state 
1  to  2  be  3  .  Then  G^n'(0  satisfy  the  following  system  of  equations 

Gi(!j(e)=d^Tlj(1'+¥TfTT4,j(^ 

4,J(«)  =  d^T^(2'  +  d?rMGi,J(£)- 


G(n+1)  m  _   V™'  .  m  +  "l["n(R'n)]  fG(n)  (E)   G(n+1)  ..... 
Gl,j  U)  -   d,(n)  +  i   V1'  +   d,(n)  +  5  [G    U)  G     U)]l,j 


.(n+1)  ...  _  A2(m-">  .  ...  +  ^min(R'"']  rG(n)  ,  .  G(n+1)  ..... 
"2,j   U)  '  d,(n)  +  i   *j(2)  +  d,(n)  +  C   [G    (5)  G     U)] 


2,j 


I, 


(3.21) 


+  dTW^TG2j1)W-  <3-22> 


+     3     r(n+D  (F\ 
+  d2(n)  +  e  Gi,j   (C) 

for  n  £  m  -  1 
where 

d,  (1)  =  A,m  +  a  ; 

d2(l)  =  x2m  +  3  ; 

d,(n)  =  A-,(m-n)  +  y-,  [min(R,n)]  +  a  ;  (3.23) 

d2(n)  =  A2(m-n)  +  y2[min(R,n)]  +  3  ; 

if  l  =  j  , 


Ui)  = 


16b. 

The  above  equations  can  in  principal  be  solved  recursively.  Finally, 

s(l.n+l)  U)   =  G(1)U)  G(2)U)  x  ^   6Cn+D(e)  .  (3.24) 

Similarly,  the  expected  first  passage  times  satisfy  the  following 
recursive  equations. 


1 


ilV1]  =  d7TTT  +  d7TTTECTi|vo  =  2] 


E[TilYo  =  2]  =  d^T  +  ^TTTECTilYo  =  1] 


E[Vilyo  =  1]  =  d^y  +  3^TE[Tn+ilvo  =  2] 


y,[min(R,n)] 
+  — 


CTHT     -{EETnlY0  =  ^  +  E^Tn+llY0  =  ^} 


T_jYn    =    2]    =   -^         +-&         EfT.  ,  I  Yft    =    1] 


n+1  >Y0  =  lj  =  d^TTTT  +  d^T  ELTn+l  'Y0 

y9[min(R,n)]    f  ^ 

+  JLT^] {ECTnJY0  =  2:  +  E:Tn+1|YQ  =  2]} 
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4.    FIRST  PASSAGE  TIME  TO  LOWER  LEVELS 

There  is  some  symmetry  in  the  process  Q,  which  we  have  not  exploited 

yet.  This  we  proceed  to  do  now.  Instead  of  the  processes  S  ,  0  <  n  <  N, 

defined  earlier,  we  now  consider  the  processes  S  ,  0  <  n  <  N. 

For  1  <  n  <  N,  S   is  the  restriction  of  the  process  Q,  observed  during 

those  intervals  of  time  spent  at  level  n,  before  the  process  Q  moves  down 

to  level  n-1  for  the  first  time.  All  S  ,  1  <  n  <  N,  are  transient  Markov 

n 

processes.  The  process  S«  is  the  restriction  of  Q  observed  at  the  lowest 
level;  it  is  an  ergodic  Markov  process.  We  denote  by  C  the  infinitesimal 
generator  of  the  process  S  ,  0  <  n  <  N. 


We  may  carry  out  for  the  processes  S   exactly  the  same  analysis  as  we 

did  for  the  processes  S  .  We  indicate  below  the  main  results  for  two 

n 

reasons. 

a.  The  first  passage  times  to  lower  levels  must  be  analysed  via  the 
processes  5  . 

b.  We  shall  be  able  to  describe  the  precise  correspondence  between  our 
analysis,  and  the  analysis  of  Neuts  for  the  infinite  quasi-birth-and- 
death  process. 


The  proof  of  the  next  lemma  is  omitted,  since  it  is  identical  to  that 
of  Lemma  2  and  Theorem  1. 


J.G 


Lemma  3 


The  matrices  C  ,  0  <  n  <  N,  are  recursively  determined  as  follows 


C  =  A^ 

Cn  =  A<nU  A^^-C^f1)  M(n+1>  ,      for  0  <  n  <  N-l. 

The  stationary  probability  vector  P^  =  (P~,  P.,  ...,  F>  )  is  determined  by 
the  equations 

^0  =  °'  (4-!) 


Pn  =  Vl  A(n"1)(-^1)J     for  Kn<N, 


N 


(4.2) 


I  P  e  =  1. 
n=0  □ 

The  (i,j)th  entry  of  the  matrix  A^n"  (~C~  )  is  equal  to 

=  (-A  n"^)  I  H'\    .  .  (4.3) 

1,1   k=l  _A(n-l)    n  K'J 

'he  factor  Ainr  /(-A;.n7  ')  is  the  probability  that,  upon  leaving  the  state 

1  ,  K        1,1 

(n-l ,i ) ,  the  process  Q  moves  to  the  state  (n,k).  The  (k,j)th  entry  of 
-C  )  is  the  expected  time  spent  by  the  process  in  state  (n,j),  starting 
from  (n,k),  before  hitting  any  state  at  level  n-l  (by  Lemma  1).  Therefore, 
:  results  from  (4.3)   that  [A^'^t-C"1)^  .  is  equal  to  (-Aj"?™)  times 
the  expected  time  spent  in  the  state  (n,j),  before  the  first  return  to  level 
n-l,  given  that  the  process  Q  starts  in  the  state  (n-l,i).  This  is  exactly 
:he  interpretation  of  Neuts'  matrix  R  for  the  infinite  quasi-birth  -and- 
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Let  t     denote  the  first  passage  time  to  level  n,  down  from  level 

m,n  3 

m+1,  for  0  <  n  <  m  <  N-l  : 


t    =  inf  {t>0  :  Xt  =  n  |  XQ  =  m+l}, 


and  let 

,(m,n)_ 


di         "ECVn   I   X0  =  m+l,Y0  =  i]    , 


d(n)   _  d(n,n) 
a.       -  a. 

for  0<n<m<N-l.  The  proof  of  the  next  theorem  is  omitted,  since  it 
is  identical  to  that  of  Theorem  2. 

Theorem  3. 

The  expected  values  of  the  first  passage  time  to  lower  levels  satisfy  the 
following  recurrence  relations. 

i'^  -  C"1  e  ,  (4.4) 

dM=  d(m,n)=  .  c-l  ,  +  ^m+D^m+l))  for  0  <  m  <  N.2         (4  5) 

—    —        m+1  —       —    '  v   ' 

dCm.n).  jCm.n+l)  +  G(m'n+1)(0)d(n)  ,  for  0  <  n  <  m  <  N-l,        (4.6) 

where   G(m'n)(0)=  n  (-Ck+1)  M(k+1),  for  0  <  n  <  m  <  N-l,  (4.7) 

k=m 


Observe  that  in  the  right-hand  side  of  (4.7)  ,  the  left  most  matrix 

"  ,) ,  not  (-  C~ ,' 
m+1'      v   n+1 


in  the  product  is  (-C"  ,),  not  (-  C  , ) 


20, 


5.  REMARKS  ABOUT  THE  COMPUTATIONAL  ALGORITHMS 

The  algorithm  A  described  in  Section  2  is  numerically  stable 
under  a  large  range  of  values  for  the  entries  of  the  generator  Q 
of  (1.1).  The  matrices  (-C~  )  have  only  non-negative  entries.  Moreover, 

since  they  measure  the  time  spent  at  level  n  only,  before  moving  to 
level  n+1,  they  usually  are  of  the  same  order  of  magnitude  for  each  n, 
even  for  large  values  of  N,  if  the  entries  of  the  matrices  A^  are  of 
the  same  order  of  magnitude  for  each  n,  and  similarly  for  the  matrices  M^n'. 
A  potential  source  of  trouble  exists  when  the  A's  are  either  very  small  or 
very  large,  compared  to  the  other  elements  of  Q.  In  that  case,  the  expected 
times  spent  at  one  level  before  moving  up  are  respectively  very  large  or 
very   small,  and  there  is  a  risk  of  encountering  overflow  or  underflow 
problems,  when  determining  the  matrices  (-C  ). 


The  steps  A3  and  A4  of  Algorithm  A  are  more  delicate.  We  start  step  A3 

with  a  vector  tt».  normalized  by  tt_n  e  =  1.  If  N  is  very  large,  it  is  likely 

that  the  vector  P..  will  be  much  smaller  than  ttm  ,  and  there  is  a  real  risk 
— N  — N 

of  running  into  overflow  problems  while  performing  step  A3.  In  order  to 
overcome  this  difficulty,  we  have  merged  the  two  steps  A3  and  A4,  and 
re-normalized  the  vectors  each  time  a  new  subvector  is  determined  (see 
Algorithm  B  in  Appendix  A). 

It  results  from  Theorem  2  that  one  may  compute  the  vectors  u}    '  '  at 
the  same  time  as  one  is  preparing  the  evaluation  of  the  stationary  probabi- 
lity distributions,  i.e.  during  step  Al  of  Algorithm  A.  In  the  numerical 
examples  presented  in  the  next  section,  we  have  determined  the  first  passage 
times  from  level  0  to  level  m,  for  1  <m  <  N.  The  corresponding  algorithm 
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Finally,  we  compare  the  numerical  efficiency  of  three  approaches. 
As  observed  by  Torrez  [12],  the  matrix  Q  is  a  band-matrix,  and  there  exist 

efficient  numerical  procedures  to  solve  eigenvectors  for  band-matrices. 

%3 
The  complexity  of  such  procedures,  for  the  matrix  Q,  is  0(N,K  ),  where 

K  =  max  {  Kn,  0  <  n  <  N)  (Wilkinson  and  Reinsh  [13] ,  p.  70).  The  crucial 

step  in  our  algorithm  resides  in  the  inversion  of  the  matrices  C  ,  0  <  n  <  N-l, 

and  the  solution  of  the  system  (2.5). 

N  3 
The  corresponding  complexity  is  0(  I  K  ).  In  Keilson  et  al.  [6],  2  N  matrices 

n=0  n 

have  to  be  inverted,  and  N  systems  have  to  be  solved.  Again,  the  complexity 

N  3 
is  0(  I  K  ).  Clearly,  the  three  approaches  have  globally  similar  numerical 

n=0  n 

efficiency;  in  order  to  distinguish  among  them,  one  would  have  to  determine  the 

coefficients  implicit  in  the  0(.)  notation.  However,  it  must  be  observed  that 

the  last  two  have  the  additional  advantage  of  offering  clear  probabilistic 

interpretations  of  the  computed  quantities. 


6.  NUMERICAL  RESULTS 

6.1  Machine  repairman  model.  We  consider  a  system  with  N=5  machines  and 
R=l  repairman.  The  system  can  be  in  K=2  environments.  The  system  remains 
in  the  environment  state  j  (j  =  1,2)  for  an  exponentially  distributed  random 
interval  of  time,  with  parameter  a-.  The  failure  rate  of  each  machine  is 
equal  to  X.  in  the  jth  environment,  with  x^   =  0.12,  and  x2  =  0.06.  The  repair 
rate  of  the  repairman  is  y  =  1  in  both  environments. 

Our  objective  is  to  measure  how  the  rate  of  changes  in  the  environment 
influences  the  system  behaviour.  In  order  to  do  so,  we  set  o^  =  6/2,  c^  =  3, 
and  chose  different  values  for  6.  Then  the  stationary  probabilities  of  being 
in  each  environment  remain  constant  and  are  given  by  yi  =  2/3  and  y^   =  1/3. 
We  expect  that  if  B  is  large,  then  the  environment  changes  rapidly,  and  the 
system  is  only  influenced  by  the  average  failure  rate  XQ  =  Yi*i+  Yox2  =  ^•** 
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On  the  other  hand,  if  6  is  small,  then  the  environment  stays  for  long 
periods  of  time  in  the  same  state,  and  this  should  affect  the  dynamic 
behaviour  of  the  system. 

We  denote  by  5^(8).  0  <  i  <  5,  the  marginal  distribution  of  the 
number  of  machines  on  repair,  for  a  given  value  of  B. 
We  furthermore  denote  by  n^X),  0  <  i  <  5,  the  probability  distribution 
for  the  classical  machine  repairman  system,  with  constant  failure  rate  x. 

In  table  I,  we  give  the  cumulative  probabilities,  corresponding  to 

-5    -2    -1 
\_{\),  for   x  =  Xq,  X1  and  x2,  and  to  £(6),  for  8  =  10  ,  10  ,  10  ,  1,  10, 

2    5 

10  ,  10  .  We  observe  that  the  distributions  of  all  the  systems  in  a  random 

invironment  are  close  to  the  distribution  for  the  system  with  a  unique, 

average  rate  Xq,  with  increasing  differences  when  the  environment  changes 

-5 
nore  slowly.  We  also  observe  that  £(10  )  =  YinUi)  +  YpH^?)*  up  t0  ^ive 

decimal  places. 

We  conclude  therefore  that  for  the  present  model,  the  random  environment 

as  little  effect  on  the  marginal  stationary  distribution  of  the  number  of 

machines  on  repair. 

In  order  to  measure  the  influence  of  the  random  environment  on  the 

lynamic  behaviour  of  the  system,  we  have  computed  the  average  time  needed  to 

reach  the  states  [n  machines  on  repair],  1  <  n  <  5,  starting  from  the  state 

machine  on  repair].  We  present  on  Figure  1  the  results  for  n=5,  which  is 

nterpretable  as  the  "time  to  complete  failure".  The  functions  are  as  follows 

1  and  2  :  f^B)  and  f2(B)»  where  f.(s)  =  E{time  to  reach  the  state  [all 

machines  on  repair]  starting  from  the  state  [0  machine  on  repair  and 

environment  j] }. 

3,  4  and  5  :  g(x,),  g(x2)  and  g(X.),  where  g(X)  =  E{time  to  reach  the  state 

[all  machines  on  repair],  starting  from  [0  machine  on  repair]}, for  a  system 

with  constant  failure  rate  i. 
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-  6  :  Yifi(B)  +  Yp^^)'  e<3ual  t0  tne  stationary  expected  time  to  complete 
failure,  starting  from  [0  machine  on  repair]. 

We  clearly  observe  several  typical  ranges  of  values  for  e.  If  6  >  10, 
the  environment  changes  so  rapidly  that  the  expected  time  to  complete  failure 
is  equal  to  the  time  for  a  system  with  unique  failure  rate  equal  to  xn. 
For  10  <  $<10,  the  expected  time  to  complete  failure  does  not  depend  on  the 
initial  environment  state.  For  e  <  10"°,  the  environment  changes  so  slowly 
that  complete  failure  is  reached  before  a  change  of  environment  occurs. 
In  fact,  it  seems  that  the  system  is  almost  completely  decomposed  in  two 

different  systems,  one  corresponding  to  each  environment,  with  \/ery   slow 

-3 
migrations  from  one  to  the  other.  Also,  we  note  that  for  &  <  10  ,  the  sta- 
tionary expected  time  to  complete  failure  looses  any  practical  significance. 


Finally,  the  average  failure  rate  \Q   yields  an  overestimation  of  the  time  to 

_3 
complete  failure  for  10  <B<1,  that  is,  for  v< 

lower,  nor  much  higher  than  the  failure  rates 


_3 
complete  failure  for  10  <&<1,  that  is,  for  values  of  $  which  are  neither  much 


6.2  Mass  search  model .  We  define  a  reference  model  with  p=15  prey  and  m 
predators.  The  system  can  be  in  K=2  environments,  with  parameters  a,  and 
a2  =  2  a,.  The  detection  and  loss  rates  are  respectively  given  by  A,  =  .001, 
X2  =  .005,  y,  =  .02  and  y2  =  -°1-  Thus,  environment  2  is  more  favorable  to  the 
predators,  since  the  detection  rate  is  higher,  and  the  loss  rate  is  lower; 
however,  environment  1  lasts  on  the  average  longer  than  environment  2. 

For  this  model,  the  rate  of  changes  in  the  environment  has  an  influence 
on  the  marginal  stationary  distribution  of  the  number  of  prey  under  surveillance 
In  Table  II,  we  give  that  distribution  for  five  different  models.  Columns  1 
and  2  correspond  to  models  with  a  unique  environment  (respectively  environments 
1  and  2). 
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lumns  3  and  4  correspond  to  the  reference  model,  respectively  with 

(slow  environment  changes),  and  a,  =  1  (fast  environment  changes). 

olumn  5  corresponds  to  a  model  with  a  unique,  average  environment  : 

2      1         2       1 

-  *2,  y  =  — j-  y,  +  -g—  vo-   N°te  tnai:  the  distribution  in  column 

has  two  modes,  respectively  corresponding  to  the  modes  in  columns  1  and  2. 

We  have  also  measured  how  the  system  responds  to  increases  in  the  effi- 

ncy  of  the  predators.  We  have  considered  two  ways  for  the  predators  to 

ime  more  efficient.  The  first  one  is  by  becoming  more  numerous,  the  second 

by  increasing  the  probability  of  detecting  a  free  prey,  or  by  decreasing  the 

robability  of  losing  a  prey  under  surveillance.  We  present  on  Figures  2  and 

3  the  values  of  t  . ,  n=5,10,15,  j=l,2,  where 
n  »j 

t  .  =  E  [time  until  n  prey  are  under  surveillance   at  time  0, 
n  ,j 

0  prey  under  surveillance,  environment  state  is  j], 

the  reference  model  with  cu  =  .0001,  and  m  equals  15  to  45.  It  clearly 

>ears  that  environment  2  is  more  favorable  for  the  predators.  The  curve  for 

presents  a  plateau  which  we  explain  as  follows.  We  denote  by  T'  .  the 

n ,  j 

until  n  prey  are  under  surveillance  given  that  initially  zero  prey  are 
•  surveillance,  for  a  system  with  a  unique  environment,  identical  to 

lment  j.  For  values  of  m  less  than  32,  say,  TJ5  ,  is  so  large  (greater 

5 

on  the  average)  that  the  reference  model  with  a,  =  .0001  switches  to 

!  before  15  prey  are  under  surveillance  (the  switch  occurs  on  the 

4 
age  after  10  units  of  time).  Once  the  model  is  in  environment  2,  it  takes 

2     3 
average  10  to  10  units  of  time  to  have  all  preys  under  surveillance. 

occurs  before  the  environment  switches  back  to  1,  and  the  total  elapsed 

4 
time  is,  on  the  average,  approximately  10  . 
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On  Figures  4  to  7,  we  present  the  values  of  t  . ,  n=5 ,  10 ,15 ,  j=l,2, 
for  models  derived  from  the  reference  model.  We  consider  a,  =  .0001, 
m  equals  to  15  and  20,  and  we  multiply  the  rates  X  .  by  vy~  ,  we  divide  the 
rates  y •  by  vy"  »  with  y  >  1»  so  that  the  probability  ratio's  "probability  of 
detecting  /  probability  of  losing"  are  uniformly  multiplied  by  y  in  each  state. 
These  figures  may  be  used  in  conjunction  with  Figures  2  and  3  to  measure  trade- 
off such  as  the  following  one.  Suppose  we  start  from  the  reference  model  with 
m=15,  and  that  we  double  the  number  of  predators  to  m=30.  This  will  entail  a 
reduction  on  t,r  ,,  which  can  be  measured  on  Figure  2.  We  can  then  measure  on 
Figure  4  that  the  probability  ratio  must  be  multiplied  by  a  factor  y  -3   in  order 
to  obtain  the  same  reduction. 

In  Table  III,  we  give  the  approximate  values  of  y  which  give  the  same  reduction 
as  doubling  the  number  of  predators,  for  m=15  and  m=20. 
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APPENDIX  A. 

Algorithm  B  :  evaluation  of  the  stationary  probability  distribution,  the 
first  and  second  moments  of  the  time  to  reach  level  m,  1  <  m  <  N,  starting 
from  level  0. 

Bl   First  passage  times  from  level  0  to  higher  levels. 
Bl.l  Initial  values. 

<  ("Cq1)  -  -(A*0))"1  ;  (Equation  (2.2)) 
u(1)  -   (-eg1)*  u^**™    ;  (3.12) 
v<«  *  Zl-t?)!LWilll'l)+    1{1)    l  (3.18) 

g(1)(o)M-c-1)a(0);g(1'1)(o)^  g(1)(o);  (3.ii) 

U^   <-  (-C^A^U*1'1)*  U™    }  (3.15) 

B1.2     Levels  2  to  N. 
for     n  =  2  to     N     do 

<  (-Cn\\)  *     V"^  »{n'l)(-C2)   a'""2')"1   ;  (2.3) 

„(n)  «.    (-C^1)(e  +  H<n-1'u(n"1)    !  (3.13) 

ud.nj  *    „(l.n-l)+  Gd.n-1)(0)  „(n)   .  (3  14) 

v(n)     -       Z(-C^!)(I  +  M<n-1Vn-1>)   u'n» 

+  (-C-^)  ^-'M"-1);  (3.19) 


^(l,n)   ^  x(l,n-l)+2u(l,n.iyn)+s(Un-l)(D)l(n).  (3>20) 

G(n)(0)  -  (-C"^)  A^"1);  (3.11) 

G(1,n)(0)  *  G^1,n"^(0)  G(n)(0);  (3.8) 

U(n)  -  (-Cn-l)(I  +  M(n'1)lj(n-1))(-C^1)  A(n_1);  (3.16) 

^l.n)  <.   u(Un'1)G(n)(.0)rf1'n"1)(I))  U(n)  }  (3.17) 

B2  Stationary  distribution. 

B2.1  Initialization  of  the  recurrence. 

(CN  *  A^l-c;!,)^-1';  (2.3) 

tt„  «-   solution  of  lMCN  =  0,Ti\.e  =  l}  (2.5) 


B2.2  Determination  of  P  ,  0  <  n  <  N. 

— n 


for  n  =  N-l  to  0  do 


(In    *   .Inxl  M^VCZ1);  (2.6) 


(re -normalize) 


N 

a   «-   lie; 
k=n  "^  ~ 


for  k  =  n  to  N  do 


{  TT^   4-   a    TT^}    } 


TABLES 


Table  I.  Cumulative  probabilities-number  of  machines  on  repair, 


Table  II.  Marginal  distribution-number  of  prey  under  surveillance 
Missing  numbers  are  less  than  5.10"  . 


Table  III.  Values  of  y  giving  the  same  reduction  as  doubling  the 
number  of  predators. 


1 
i 

0 

1 

2 

3 

4 

5 

n(x0) 

.56395 

.84593 

.95872 

.99256 

.99932 

C(105) 

.56395 

.84593 

.95872 

.99256 

.99932 

c(io2) 

.56398 

.84589 

.95870 

.99255 

.99932 

5(10) 

.56421 

.84561 

.95850 

.99249 

.99931 

5(1) 

.56579 

.84397 

.95711 

.99202 

.99925 

^io"1) 

.56910 

.84151 

.95424 

.99091 

.99908 

5(io"2) 

.57033 

.84078 

.95321 

.99047 

.99900 

dio"5) 

.57050 

.84068 

.95306 

.99040 

.99899 

n(xx) 

.49516 

.79226 

.93486 

.98620 

.99852 

n(x2) 

.72118 

.93754 

.98946 

.99881 

.99993 

Table  I 


Cumulative  probabilities-number  of  machines  on  repair 


n 

I 

II 

III 

IV 

V 

0 

.00073 

.00048 

1 

.00820 

.00543 

.00005 

.00003 

2 

.04019 

.02662 

.00061 

.00043 

3 

.11319 

.07504 

.00428 

.00342 

4 

.20375 

.00002 

.13524 

.01959 

.01724 

5 

.24653 

.00026 

.16400 

.06175 

.05835 

6 

.20544 

.00220 

.13776 

.13712 

.13601 

7 

.11886 

.01270 

.08407 

.21606 

.22012 

8 

.04755 

.05082 

.04948 

.24021 

.24628 

9 

.01294 

.13835 

.05546 

.18523 

.18752 

10 

.00233 

.24903 

.08481 

.09610 

.09441 

11 

.00026 

.28299 

.09427 

.03198 

.03001 

12 

.00002 

.18866 

.06255 

.00633 

.00560 

13 

.06531 

.02160 

.00066 

.00054 

14 

.00933 

.00308 

.00003 

.00002 

15 

.00031 

.00010 

Table  II 


Limiting  distribution-number  of  prey  under  surveillance 
Missing  numbers  are  less  than  5.10 


m 

*5.1 

t10,l 

t15,l 

%*2 

t10,2 

t15,2 

15 
20 

4.1 
4.1 

4.4 
4.0 

7.8 
3.9 

4.9 

4.5 

7.1 

5.6 

>  15 
5.4 

Table  III 


Values  of  y     giving  the  same  reduction  as  doubling 
the  number  of  predators 


FIGURES 


Fig.l.   Expected  time  to  reach  the  state  [all  machines  on  repair], 
starting  from  [no  machine  on  repair]. 

Fig. 2.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  1], 
varying  number  of  predators. 

Fig. 3.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  2], 
varying  number  of  predators. 

Fig. 4.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  1], 
15  predators,  varying  probability  ratio. 

Fig. 5.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  2], 
15  predators,  varying  probability  ratio. 

Fig. 6.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  1], 
20  predators,  varying  probability  ratio. 

Fig. 7.   Expected  time  to  reach  the  states  [n  prey  under  observation], 
starting  from  [0  prey  under  observation,  environment  2], 
20  predators,  varying  probability  ratio. 
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Fig. 2.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 

from  [0  prey  under  observation,  environment  1],  varying  number  of  predators 


Fig. 3.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 

from  [0  prey  under  observation,  environment  2],  varying  number  of  predators 
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Fig. 4.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 
from  [0  prey  under  observation,  environment  1],  15  predators,  varying 
probability  ratio. 
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Fig. 5.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 
from  [0  prey  under  observation,  environment  2],  15  predators,  varying 
probability  ratio. 
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Fig. 6.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 
from  [0  prey  under  observation,  environment  1],  20  predators,  varying 
probability  ratio. 
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Fig. 7.  Expected  time  to  reach  the  states  [n  prey  under  observation],  starting 
from  [0  prey  under  observation,  environment  2],  20  predators,  varying 
probability  ratio. 
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